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The motion of point vortices with periodic boundary conditions was studied by using Weier- 
strass zeta functions. The scattering and recoupling of a vortex pair by a third vortex becomes 
remarkable when the vortex density is large. The clustering of vortices with various initial con- 
ditions is quantitated by the L function used in the point process theory in spatial ecology. It 
is shown that clustering persists if the initial distribution is clustered like an infinite row or a 
checkered pattern. 
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The statistical approach to the problem of assemblies 
of point vortices (PVs) dates back to Onsager (1949). 
A state of negative temperature is considered to be re- 
lated to the clustering of vortices rotating in the same 
direction and the inverse energy cascade predicted in 
the two-dimensional Navier-Stokes (2D NS) turbulence. 
In many numerical simulations, PVs are bounded in a 
circular wall since a velocity field due to a PV can be 
computed by including a single mirror image. Although 
the axisymmctry with respect to the origin is conserved, 
spatial homogeneity is not guaranteed in such a circular 
system. The numerical difficulty in the simulation of vor- 
tices in a box is that there emerges an infinite sequence 
of virtual images. Our objective is to study the turbu- 
lent motions and clustering of many PVs in a periodic 
box using Weierstrass zeta functions. 

Let us start by representing the 2D NS equation in 
terms of a complex position z — x+iy, velocity q = u—iv, 
pressure p, and kinematic viscosity v as 



qt + qqz 



-2p z + ivq z 



(1) 



Here, q denotes the complex conjugate of q and we use 
the relations d x = d z + dz, d y = i(d z — dz), u = (q + 
q)/2, v = i(q — q)/2, and A = Ad zz . The incompressible 
condition yields V • v — q z + q z = 0. The vorticity u> = 
v x — u y can be expressed by q as u = 2iq s . 

If the flow is irrotational lu — 0, then q z = 0, q depends 
on only z (and t), and the theory of conformal mapping 
can be applied. Equations for the vorticity and pressure 
are 

Qrt + QQzm + Mzz = ^Qzzz, (2) 

and p zz = -(q z q z + q z )/2, respectively. 

According to Tkachenko, 1,2 the velocity field due to 
a single PV at the origin with periodic boundary condi- 
tions (BCs) is equivalent to that due to the PVs on the 
lattice z mn = 2mwi + 2nuj2 , where the complex numbers 
uj\, 0J2 are the half periods of the lattice and to, n are ar- 
bitrary integers. The ratio of the two periods r = u>i/u>2 
can be restricted in the region 



Imr > 0, |Rer| < 1/2, \t\ > 1. 



(3) 



We investigate the case of square periodic BCs, which 
are usually employed in the numerical studies of two- 
dimensional turbulence by selecting r = i. However, we 
can deal with an arbitrary periodic parallelogram by con- 
sidering various values of r that satisfy (3). 

The velocity field due to a PV of strength k = 2ir is 
given by the Weierstrass zeta function £(2; 0J2) along 
with a rigid rotation term as follows: 



q — iQ(z) — iilz = w(z). 



(4) 



Since the vortex lattice undergoes rigid rotation with 
an angular velocity f2 = 7r/[4Im((DitJ2)], the second 
term in Eq. (4) is necessary in order to cancel the ve- 
locity circulation on the boundary. The vortex density 
n = l/[4Im((DiW2)], ^, and the vortex strength k are re- 
lated as kit, — 2fL If the length of the side of the square 
is unity, then u>i = 1/2, u>2 = */2, f2 = 7r, and k = 2n. 

The equation for the streamline tf> — const., where if) 
is the streamfunction, is equivalent to dx/ip y — —dy/ipx. 
Using the relations u = ip y and v = —ip x , ^ is expressed 
as ip — Judy + fix). The sigma and zeta functions 
of Weierstrass are related as £(z) = a' (z)/a(z). This re- 
lation is consistent with the asymptotic forms £ ~ 1/z 
and a ~ z when z ~ 0. Using the above results, ip for a 
single vortex lattice centered at the origin is given by 



ip = —Rein a(z) 



/2- 



(5) 
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There is a minimum value of ip in the periodic case, in 
contrast with the unbounded plane in which there are no 
limits on ip . 

For simplicity, we consider an assembly of PVs with 
Ki = 2n(ii, fii = 1, for i — 1, • • • , N% and \ii = —1 for 
i = Ni + 1, • • • , N(= Ni + N 2 ). Therefore, ip for N PVs 
located at zt is given by 

V> = S^ lA/l {-Re[lna(z - z t )} + Q\z - Zl \ 2 /2}. (6) 

Using Eq. (4), the equation of motion of PVs with square 
periodic BCs can be expressed as 3 

Zi = EjjLiiijwizi - zj). (7) 

The equation can be rewritten in the Hamiltonian form 
as fiidzi/dt = dH/dzi, where the Hamiltonian H can be 
expressed as H 



EiLi Mi 



EN sr^N -l 
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Fig. 1. Trajectories of three (a) and four PVs (b). PVs 1, 2, 3, 
and 4 are denoted by solid, dashed, dotted, and dotted-dashed 
curves, respectively. 

hij = -Rc[ln a(z l - Zj)] + Q\zi - Zj\ 2 /2. The Hamilto- 
nian, which is given by the total kinetic energy minus the 
self-induced kinetic energy of PVs, can be interpreted as 
the sum of a kinetic energy due to interactions between 
the pairs of PVs. 

If Mathematica is used, we can compute the Weier- 
strass zeta function as we use the sinusoidal function in 
a Fortran code. The system (4,7) is solved numerically by 
the NDSolve command of Mathematica 5.2 installed in a 
PC with an AMD Athlon 64x2 3800 CPU, 2 GB memory 
and Windows XP OS. The computation is realistic since 
the CPU time in such a PC environment ranges from two 
to five days for 100 PVs and 10 eddy turnover times. 

If the PVs lie in an unbounded domain, the system has 
four integrals: 4 the Hamiltonian H u = — ^ p-iPj In \z% — 
Zj\, two components of the linear impulse I = (I x ,I y ), 
I x = ^/UjRe[zi], I y = ^/Xjlmfzj], and the angular im- 
pulse A — Y^, Pi\zi\ 2 - Since the system in a periodic box 
has no circular symmetry, A is no more constant; how- 
ever, H,I X , and I y are conserved. Since there are three 
conserved quantities, the system of three PVs is inte- 
grable, while the four PVs exhibit chaos. 

Examples of the trajectories of three PVs with fa = 
(2,2,-1) and an initial condition (21,22,23) = (0,0.5, 
0.25 + i\/3/4) located at the vertices of an equilateral 
triangle and that of four PVs with fa = (2,2,-1,-1) 
and (zi,z 2 ,z 3 ,2 4 ) = (0, 0.5, iV3/4, 0.5 + i\/3/4) at t = 
located at the vertices of a rectangle are shown in Figures 
la and lb, respectively. The first case leads to a collapse 
if the PVs are in an unbounded plane. 

To analyze the spatial distribution of many PVs, we 
introduce the L function used in the point process theory 
in spatial ecology. 5 ' 6 Let afj = (xi,yi) be the position 
of N points randomly distributed in an area S. The K 
function is defined by 

K(r) = {\N)-^ti^U^ r ~ l*i " ^1). ( 8 ) 
where A = N/ S is the number density of the points and 
6{x) is the step function. An extra function added in (8), 
in order to modify the edge effect, is unnecessary in the 
present periodic case. 

If the distribution of points obeys completely spatially 
randomness (CSR), which is synonymous with a homo- 
geneous Poisson process, K (r) becomes the area of the 
circle with radius r, i.e., K(r) = irr 2 . Then, it is conve- 



nient to introduce the L function as 

L(r) = y/K(r)/ir-r. (9) 

CSR yields L = 0. For clustering, i.e. points staying close 
to the other points, we have L > 0. If the points tend to 
be at a distance from each other, L < 0. Whether L{r) 
is positive or negative can depend on r. For instance, 
a checkered pattern (Ichimatsu moyo in Japanese) gives 
L > for small r, but L < for large r. 

According to Novikov (1976), 7 for an unbounded 
plane, we have a relation between the distance rji of two 
PVs of strength Kj and m and the energy spectrum E(k) 
as 

E{k) = (47rfc) _1 [£jK* + 2X j<l K j KiJ {kr j i)}. (10) 

If Kj — k for all j, we have 

E(k) = K 2 (47rfc)- 1 [V + 2S J<J J (fcr j7 )]. (11) 

Using the number density p(r) at the distance r be- 
tween two PVs, we have, in the continuous limit, 
E(k) = E^k) + E 2 (k), where E^k) = k 2 V(4tt/c)- 1 
and E 2 (k) = k 2 (2tt/c)- 1 J °° J (kr)p(r)dr. Ei(k) and 
E 2 (k) correspond to the self-energy of each vortex and 
the interaction energy between two PVs, similar to hij. 
The relation between the K(r) and p(r) is K(r) = 
PN'^ 1 J Q r p(r')dr' , where p(r) is normalized by the up- 
per limit / of r and the total number N' = N(N - l)/2 
of pairs of PVs. 

If p(r) = Cr a , we can integrate E 2 {k) into E 2 (k) = 
n 2 Ck- a - 2 r((a + l)/2) ^((l-a)^)]- 1 for -1< a < 
1/2. The value a = —1/3 gives the Kolmogorov spectrum 
m a 3D turbulence. 

Formally, the CSR value a = 1 in two dimensions 
yields the fc~ 3 spectrum, although the integral does not 
converge. In order to avoid divergence, an exponential 
decay in p(r) may be introduced. Otherwise, we may set 
the upper limit I in the range of integration in E 2 {k). Us- 
ing the normalization r = lr' and the integral expressed 
by a regularized hypergeometric function as J Q Jo(fcr) 
r a dr = r((l + a)/2) 1 F 2 ((l+a)/2; {1, (3+a)/2}; -fc 2 /4) 
[2r((3 + a)/2)]-\ for a > -1, E 2 (k) is written as 
E 2 (k) = n 2 Cr +1 T((l + a)/2) iF 2 ((l + a)/2; {1, (3 + 
a)/2} ;-(fc0 2 /4) [47rfcr((3 + a)/^]" 1 . 

If we use the CSR distribution p(r) = nN'l 2 r, we 
have Jo(kr)rdr = J\{k)/k. Then, E 2 (k) is given by 

E 2 (k) = K 2 NT 1 k- 2 J 1 (kl)/2. (12) 

For large kl, we have the asymptotic form 

E 2 (k) ~ y/l/2TVK 2 N'l(kl)- 5/2 cos(kl - 3tt/4). (13) 

The value of E 2 (k) is oscillatory as k increases with the 
amplitude decaying as fc -5 / 2 . Evidently, the total energy 
spectrum docs not become negative because E\{k) 3> 
E 2 (k) for large k. 

For numerical studies, we first consider an assembly of 
PVs having the same positive strength k (= 2ir). The 
following four typical cases are considered: Case (I) an 
infinite row thatis a discrete model of the vortex sheet, 
Case (II) PVs located randomly in checkered patterns, 
Case (III) PVs located randomly in the 10 x 10 sub- 
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Fig. 2. The L functions, (a) corresponds to (I), (b) to (II), (c) to Fig. 3. The L functions, (a) corresponds to (V), (b) to (VI), and 



(III), and (d) to (IV) 



(c) to (VII). (d) shows the distribution of PVs for (VI) at t = 0.1. 
White (black) circles denote positive (negative) PVs. 



squares, and Case (IV) CSR in the unit square. Here, the 
word CSR implies that the PVs are distributed by using 
random numbers generated by a single run. The initial 
conditions, the number TV of PVs, the final time tf, and 
the values of three conserved quantities are summarized 
in Table 1. The relative precision of H in the numerical 
simulation is confirmed to be less than 10~ 6 ~ 10~ 5 up 
to t = tf. Figure 2 shows the L function computed by 
the initial and final distributions of PVs. 

For Case (I), the PVs are initially located on the x- 
axis as Zj(0) = j/N+ esm2nj /N, j — 1, • • • ,7V, where 
e = 0.05. If e is fixed and N is increased, the pairing of 
two adjacent PVs becomes more conspicuous than the 
winding of the sheet due to the Kelvin-Helmholtz insta- 
bility. The growth rate a = mrp(l — p)/a 2 of the pairing 
instability in an unbounded plane attains a maximum at 
the wavcnumbcr p = 1/2, where a is the distance be- 
tween two adjacent PVs. 4 Modulus 1 is suitably consid- 
ered so that the PVs are plotted in the selected square. 
Of course, the PVs wander chaotically from one square 
to another. 

For e = 0, we have p(r) = 2N'T 2 , K(r) = 2r 
and the average Hamiltonian (h) = 2^ li> j hij[N(N — 

I)]- 1 = SEfii" 1 i>(j/N)+iP(l/2)/2][N(N- I)]- 1 ~ 

1 /2 

2 L dxip(x). Since ip(x) <~ — ln|x| for x ~ 0, (h) = 
J s dxdyip(x,y)p(r) yields a finite value for p = Cr a , 
r ~ 0, and a > —2. Since a < 1 corresponds to clus- 
tering, clustering with a < — 2 cannot occur from the 
initial condition with a finite (h). 

Second, we consider Case (II) where the initial PVs are 
located randomly in eight segments revealing a checkered 
pattern. L(r) > for < r < 0.15 implies that the 
PVs are clustered, while L(r) < for 0.15 < r < 0.4 
means that they are uniformly spaced at larger scales. We 
observe that this tendency remains at t — tf, although 
it becomes somewhat weak and an additional oscillatory 
behavior is observed. 



Third, Case (III) initially has uniformly spaced PVs. 
L(r) has a negative value and a minimum at r s=s 0.06. We 
observe that L{r) is almost the same during t € [0, t/]. 

Fourth, we examine the completely spatially random 
distribution of the PVs at t = in Case (IV) . We observe 
slightly negative values of L(r) at t = and tf, but we 
do not see any significant differences between the initial 
and final distributions. 

Next, we consider the system with both positive and 
negative PVs of the same strength and the same num- 
bers Ni — N 2 = N/2. The following three typical cases 
are examined: Case (V) the Karman vortex street, Case 
(VI) positive and negative PVs located alternately in 
checkered segments (16 subsquares), and Case (VII) com- 
pletely spatially random distribution. The initial condi- 
tions, the number of PVs, tf, and the values of three 
conserved quantities are summarized in Table 2. 

First, we consider Case (V) given by the follow- 
ing expression: (pj,zj(0)) = {—l,j/N\ + e(R + iR)) 
for j = 1,... ,N X and (^,«,-(0)) - {l,j/N 2 + 1/N+ 
ih + e(R + iR)) for j = iVi + 1, • • • , N, where e = 10~ 4 , 
R denotes random numbers, and the distance h between 
two rows is taken as 1/N. In this case, the negative PVs 
are first numbered. At approximately t = 0.001, we ob- 
serve that the two rows begin to break with the pairing 
instability. 4 

For the two types of PVs, we introduce Ki m (r) for 
{I, m) = (+, +), (-, -), and (+, -) as 

K lm (r) = (XN\)~ 1 T,iT,j6(r - \x t - xj\), (14) 

where the sum is for i, j = 1, • • • ,N\ except for j = i 
if (l,m) = (+,+)• The case (l,m) = (— ,— ) is similar. 
The sum for (I, m) = (+, — ) is both for i = 1, • • • , Ni, 
j = JVi + 1, • • • , AT and i = JVi + 1, • • • , N, j = 1, • • • , JVi. 
We define Li, n (r) from Ki m (r) similar to Eq. (9). We 
observe strong clustering for (/, m) = (+, — ) at t = 0.01 
corresponding to the pairing instability. 
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Second, for Case (VI), we observe the asymmetry, i.e., 
L++(r) > L (r) for all r at t = 0.1. The initial clus- 
tering is incidentally stronger for positive PVs than the 
negative ones. This fact can be regarded as clustering to 
a single vortex at the largest scale. There is also a signif- 
icant void where positive (or negative) PVs do not exist 
at t = tf. We can also consider the initial fractal distri- 
bution like a Sierpinski's gasket, which shows clustering 
at large scales. 

Finally, we consider the CSR distributions of two types 
of PVs as that in Case (VII). A remarkable feature in 
this turbulent situation is that there are several pairs 
of positive and negative vortices moving linearly at a 
velocity of k/4/wt, where 2h is the distance between two 
PVs. Since the pair is surrounded by a number of other 
isolated PVs, however, the moving direction is bent by a 
third vortex when they cross each other. Moreover, if the 
collision is nearly head-on, a vortex of the pair with a sign 
opposite to the third target vortex replaces its partner 
with the latter and then continues to move linearly again. 
An exact analysis of such scattering of three PVs in an 
unbounded domain was already given by Aref (1979). 8 

The examples of scattering and recoupling of three 
PVs in a periodic box are given by an initial location 
(zi, Z2, Z3) = (L + iL, (L + d + h)i, (L + d — h)i) with 
L = 1/2. A pair of vortices 2 and 3 is initially approach- 
ing vortex 1. Figure 4a shows their trajectories when 
h = 0.02 and d = ±0.02, ±0.01, 0.04, 0.08, 0.16, and 0.32. 
The final time tf is 0.04 except for tf = 0.1 for d = 0.16, 
tf = 0.05 for d = 0.32, and t f = 0.06 for d = -0.01. 
Recoupling is observed for d = 0.01,0.02,0.04,0.08, and 
0.16. The dependence of a 7r-normalized scattering angle 
S(j)/ir measured by the moving direction of z 3 (t = 0.04) 
for h = 0.02 in a periodic box is shown in Figure 4b. The 
recoupling of PVs in an unbounded plane with L — > oo 8 
occurs if < d/h < 9. On the other hand, the present 
simulation in a periodic box with h = 0.02 shows a shift 
of the range d for recoupling as —0.5 < d/h < 8.5, al- 
though this range may vary as h is changed in the case of 
periodic BCs. In a GIF animation, successive scattering 
and recoupling, similar to the chaos in a billiard system, 
are clearly observed. The existence of such vortex pairs 
may play a crucial role in stirring assemblies of PVs. 

Since the average distance of randomly located PVs 
is I ~ N~ x l 2 and the strength 2ir is fixed, the typical 
velocity and eddy turnover time are v ~ N 1 / 2 and t e ~ 
1/N, respectively. Denoting the smallest distance of the 
vortex pair by al, its velocity is V <~ 1/al. Because of the 
recoupling condition, the cross section of the scattering 
is approximately a c <~ d <~ Wal and the area swept by 
the pair during t e is S <~ dVt e ~ 10/N. Therefore, the 
condition for the scattering to occur in t e is N <~ 10 since 
S <~ 1, the size of the square. If N = 100, t e ~ 0.01 and 
one pair will be scattered approximately 10 times in a 
numerical simulation in the time interval < t < 0.1 ~ 

The Li m (r) function becomes slightly positive for 
(l,m) = (+,—), which also indicates that the pairs of 
positive and negative PVs survive until t = 0.1. How- 



ever, the absolute values are much smaller than the ini- 
tially clustered cases. To clearly observe the spontaneous 



Table I. Single type of PVs 



Case 


N 


*/ 


H 


Ix 


h 


I 


100 


0.01 


17948 


50.5 





II 


96 


0.1 


11973 


47.872 


47.923 


III 


100 


0.1 


12796 


49.968 


50.020 


IV 


100 


0.1 


12886 


49.687 


50.206 



Table II. Two types of PVs 



Case 


iVi, 2 


*/ 


H 


Ix 


h 


V 


50 


0.01 


-415.31 


-0.50019 


-0.49976 


VI 


48 


0.1 


369.60 


-12.176 


-11.081 


VII 


50 


0.1 


-178.31 


-1.9433 


2.3455 




Fig. 4. (a) Trajectories of three scattering and recoupling PVs for 
various values of h. Solid, dashed, and dotted curves denote zi, 
Z2, and 23, respectively, (b) The Tr-normalizcd scattering angle 
S<f>/n versus d for h = 0.02. 

clustering, longer simulations may be required. We also 
investigated the probability distribution function of ve- 
locity circulation, which is studied in Umeki (1993) 9 for a 
3D turbulence. 10 A similar approach in the point process 
theory is called the Quadrat method. 5 

In summary, a method to simulate the motions of PVs 
with periodic BCs is described. Several numerical exam- 
ples are illustrated and the clustering of PVs with differ- 
ent conditions is examined by the L function. 
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